Statistical Mechanical Calculation of Anisotropic Step Stiffness of a Two-Dimensional 
Hexagonal Lattice Gas Model with Next-Nearest-Neighbor Interactions: Application 

to Si(lll) Surface 
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We study a two-dimensional honeycomb lattice gas model 
with both nearest- and next-nearest-neighbor interactions in 
a staggered field, which describes the surface of stoichiomet- 
rically binary crystal. We calculate anisotropic step tension, 
step stiffness, and equilibrium island shape, by an extended 
random walk method. We apply the results to Si(lll) 7x7 
reconstructed surface and high-temperature Si(lll) lxl sur- 
face. We also calculate inter-step interaction coefficient. 

PACS numbers: 68.35.Md, 68.35.-p, 50.50.+q, 64.60.-i 



I. INTRODUCTION 

Recent developments of the experimental techniques 
such as STM (scanning-tunneling microscopy) jl| , LEEM 
(low-energy electron microscopy) || and REM (reflection 
electron microscopy) || make it possible to observe a step 
on a crystal surface in the wide range of the length scales. 
However, the connection among quantities measured in 
different scales has not been clarified yet. 

In Ref. 0, for two-dimensional (2D) square-lattice 
Ising model with both nearest- and next-nearest-neighbor 
(nn and nnn) interactions, we calculated anisotropic 
interface tension and interface stiffness by the imagi- 
nary path-weight (IPW) method which is an extended 
Feynman-Vdovichenko's random walk method |^-||. In 
the method, the overhang structure in a step is taken 
into account, which leads to high accuracy of the results 
in a wide range of temperature. 

We applied the results to Si(OOl) surface based on the 
microscopic kink energy obtained by the Swartzentruber 
et al. H . The Ising result gave a satisfactory explanation 
for experimentally measured step tension 7, step stiffness 
7 and equilibrium island shape obtained by Bartel t et 
al. @ on Si(OOl) surface. 

In the present paper, we consider the honeycomb lat- 
tice Ising system in a staggered field, with both nn and 
nnn interactions, to calculate interface tension, interface 
stiffness, island shape and the coefficient of step inter- 
action by the IPW method. We aim at applying the 
results to Si(lll) 7x7 reconstructed surfaces and high- 
temperature Si(lll) lxl surface. 



II. MODEL HAMILTONIAN 

We consider a honeycomb lattice with 2N sites. We 
decompose the lattice into two triangular sublattices des- 
ignated by A and B. On the A-sublattice, we define the 
occupation variable Cm which takes 1 (present) or (ab- 
sent) at the site i. Similarly, we define C B j for the B- 
sublattice. 

The lattice gas Hamiltonian 7Ylg is then written as 
"Hlg = — 4Ji ^ [CajCbj — 2^ Ai + ^Bj)] 

<i,j> 

— 4Ja2 ^ [CajCa 

<i,j> 

— 4JB2 22 [^BiCb. 

N 

" e A @Ai — £B CbU (2-1) 



<i,3> 



-(C Al + C A j)} 



^(CBi + C Bj )] 



N 



where AJ\ is the bond energy between A-atom and B- 
atom of the nn sites, 4Ja2, and 4Jb2 are the bond ener- 
gies between nnn atoms. ca(£b) is the "surface chemical 
potential" of A-atom (B-atom) . We consider the simplest 
case where £a and £b are given by 



£A = Ma, gas (-Pa , Pb , T) - ha, S ux{(T), 

£B = MB,gas(-PA, ft, T) - ^B.surfCP)- 



(2.2) 



In the above, /JA,gas(PA, Pb, T) is the chemical poten- 
tial of A-atom in the gas, Pa is the partial pressure 
of A-atom in the gas phase (similarly for MB,soiid(P), 
MB,gas(pA,-PB,P) and P B ); MA, SU rf(P) and ^B,surf(P) are 
expressed as 



MA.surf = A-E(T) + AWidCO, 

MB.surf = - AE(T) + HsoUd(T), 



(2.3) 



where AE(T) has been introduced as the difference 
from the chemical potential of atoms in the bulk solid 

A*solid(P)- 

Let us consider the bulk phase-coexistence state of the 
stoichiometrically binary system. Total chemical poten- 
tial of the system has to be unchanged under removal of 
one pair of AB atoms from crystal into vapor and vice- 
versa. Hence, as the coexistence condition, we have 
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= 2 Mso i id (T). (2.4) 



Combining (2.2)- (2.4), we obtain 
£b = -£A- 



(2.5) 



This condition also means that, in the lattice gas Hamil- 
tonian (2.1), the total energy of the all-occupied state is 
the same at that of the all-empty state. 

Let us introduce the Ising spin variables {a A i} and 
{<7 B i} as 



< T \, = - ( Cm — - j , ctbi = 2 ( - — Csi 



(2.6) 



Substituting (2.6) into the Hamiltonian (2.1) together 
with (2.5), we have the Ising AF Hamiltonian Ti: 

H = H LG + Nz 6 Ji + Nz 3 (J A2 + J B2 )/2, 

= Jl ^ PAiCBj ~ JA2 ^ &Ai&Ai 



<i,J> 



<i,3> 
N 



N 



— JB2 ^ PBiVBj — H / ] CTAi ~ H ^ OBi 
<i,j> i—1 i — 1 

+Nz 6 J 1 + Nz 3 {J A2 + J B2 )/2, 
H = e A /2 (2.7) 

where zg = 3 and z$ = 6 are the coordination numbers 
of honeycomb lattice and triangular lattice, respectively. 

III. IMAGINARY PATH- WEIGHT METHOD 

We calculate interface quantities by a random walk 
method with imaginary path- weight (IP W) |^-||]. We re- 
gard an interface with zigzag configuration as a trace of 
2D free random walk. 

Consider an interface which connects site O and P. 
We denote the distance between site O and P by R (Fig. 
|l|) . The interface of the two-dimensional Ising model is 
made by fixing the boundary spins as depicted in Fig. [I]. 
Let us denote the partition function the Ising model with 
and without interface by Z~^~ {9) and Z~^ + , respectively, 
where 9 is the slant angle of an interface relative to a lat- 
tice axis. |IH , |i2[ ] The interface tension j(6, T) is defined 
as 



1 (6,T) = -k B T lim i In 

fl-»oo R 



7+' 



(3.1) 



where Z~^~ (9) / Z~^ + is regarded as the interface partition 
function Q . 

We apply Vdovichenko's method |j| to deal with 
the low-temperature diagrammatic expansion of Zi^(0) 
and Z^ + . The method, which originally treated the 
high-temperature expansion of the partition function, 



also works for low-temperature expansion to evaluate 
weighted sum over all possible domain-wall configura- 
tions. The essential point in the Vdovichenko's method 
is introduction of the imaginary factor e 4 */ 2 at each turn 
(with angle <p) of the random walks of the domain wall. 
With this simple recipe, the problem reduces to a free 
random walk problem on a lattice. 

We see that Z^ + equals weighted sum over possible 
configurations of closed domain walls; and Z^~(9) equals 
weighted sum over possible configurations of closed do- 
main walls plus a single "open" domain wall travers- 
ing the lattice. In evaluating Z^~(9) by Vdovichenko's 
method, the free random walk nature allows us to "de- 
couple" the open domain wall from closed domain walls. 
H(?J Therefore, in the limit of R — > oo, the interface par- 
tition function is equivalent to the "edge-to-edge" lattice 
Green's function of the free random walk with running 
on the dual lattice. JjJ Thus, in the limit of R — > oo, the 
interface partition function Q is written as [jjj 



Q = cxp(-7(0)i?/fc B T) 



gifc-R 

dkxdky D(k)> 



where the D-function is defined as 



D(k) = det[l - A(k)]. 



(3.2) 



(3.3) 



Here, A(k) is the Fourier component of the connectivity 
matrix A(r) which characterizes the random walk. 

The above-described imaginary path-weight random 
walk method to calculate 7(0, T) is exact only for solv- 
able cases. However, we have verified that the method 
works fairly well also for non-solvable cases PJ^, |l3| , |l4| . 

After evaluating the integral by the pure imaginary 
saddle point ui = (u> x ,ojy), we obtain a set of equations 
as 



Diius) = 0, 

dD(ju) ,dD(iw) 

-/ „ =tan( 



dUJy 



and 



j{9,T) = k B T(io x cosfl 



(3.4) 



(3.5) 



From the thermodynamical theory on equilibrium crystal 
shape(island shape) |l6| , we have, 



= Xy/k B T, 



Xx/k B T, 



(3.6) 



where A is the Lagrange multiplier associated with the 
volume-fixing constraint in the Wulff construction, and 
x and y are the Cartesian coordinates describing the 2D 
island sha pe. Th us, we obta in t he island shape directly 
from (3.4) with (3.5). Eq. (3J3) gives relation between 
the interface orientation angle 9 and the point (x, y) on 
the island shape. 
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The interface stiffness, which we denote by j(6), is 
given by |J| 



1(0) = 1(0) 



dff 2 



k B T^D 2 x +D 2 y -[~D xx sin 
+D xy sin 29 - D yy cos 2 9} ~ 1 , 



(3.7) 



where 



D _ 3D _ 3D 

x 8uj x ' " diu x ' 

D D 

xx duv m 8lu?: 



n ^2- 

OLU X Uy 



(3.8) 



The one-dimensional interface of the lattice gas corre- 
sponds to a step on the vicinal surface. For the vicinal 
surfaces, the surface free energy per projected area, which 



we denote by f(p), is written as [17 



19 



f(p) = f(0) + 7(9) P + B(9)p 3 

p = — tan<^>, 
ah 



(3.9) 



where p is the step density, <j> is the tilted angle of the 
vicinal surface, and et^ is the step height, and B(9) is the 
inter-step interaction coefficient. We have |20| , pl|| 

B(0) = 4^( 90 ), 



(3.10) 



A(ffo) = g ( 1 




where go is the coupling constant of the long range inter- 
action between the steps of the form go / r 2 (r s is the step 
separation distance). Note that, in the limit of go — > 0, 
the factor X(g) approaches unity, leading to |l§| , p| 



B{9) 



tt 2 (fc B r) 2 

6 y(6) 



Hence, the stiffness (3.7) can be utilized in determining 
the inter-step interaction coefficient B(9). 



IV. THE CONNECTIVITY MATRIX AND THE 
D FUNCTION 

A. The honeycomb lattice gas model with next 
nearest neighbor interaction 

We apply the IPW method to calculate interface quan- 
tities of the nnn Ising model on the honeycomb lattice 



described by Hamiltonian (2.1,2.7) (Fig. ||). The Fourier 
components of the connectivity matrix {A mn ) are 



Ma 


= 


exp 


ik x )W/W H , 




A 2 ,i 


= 


exp 


ik x )W/W H r pi 




A 3 ,i 


= 


exp 


ik x )W/W H r p r p W a , 




A$,i 


= 


exp 


ik x )W/W H r m r m W b , 




Ae,i 


= 


exp 


ik x )W/W H r m , 




A h2 


= 


exp 


ik x /2 + \k y c y )W ■ W H r m , 




A 2 ,2 


= 


exp 


ik x /2 + ik y c y )W-W H , 




Az,2 


= 


exp 


ik x /2 + \k y Cy)W ■ W H r pi 




A 4 ,2 


= 


exp 


\k x /2 + \k y Cy)W ■ W H r p r p 


w b , 


A S ,2 


= 


exp 


\k x /2 + \k y Cy)W ■ W H r m r 


mW a , 


Ai, 3 


= 


exp 


-ik x /2 + ik y c y )W/W H r m 


r m W bl 


A 2 . 3 


= 


exp 


-ik x /2 + ik y Cy)W/W H r m 




A 3 ,3 


= 


exp 


-ik x /2 + iky C y)W/W H , 




Ai, 3 


= 


exp 


-ik x /2 + ik y Cy)W/W H r p , 




A 5 ,3 


= 


exp 


-ik x /2 + ik y Cy)W/W H r p r 


P w a , 


A 2A 


= 


exp 


-ik x )W ■ W H r m r m W a , 




A 3 ,A 


= 


exp 


-ik x )W ■ W H r mi 




A 4A 


= 


exp 


-ik x )W-W H , 




A 5A 


= 


exp 


-\k x )W ■ W H r p , 




A 6A 


= 


exp 


-ik x )W ■ W H r p r p W b , 






= 


exp 


-ik x /2 - ik y Cy)W/W H r p r 


pW a , 


A 3 , 5 


= 


exp 


-ik x /2 - \k y Cy)W/W H r m 


r m w b , 


Ai, 5 


= 


exp 


-ik x /2 - \k y Cy)W/W H r m 




A 5 , 5 


= 


exp 


-ik x /2-ikyCy)W/W H , 




A 6 ,5 


= 


exp 


-ikx/2 - ik y Cy)W/W H r p , 




Al, 6 




exp 


ik x /2 - \k y c y )W ■ W H r p , 




A 2 ,6 




exp 


ik x /2 - \k y c y )W ■ W H r p r p 


w b . 


A 4 ,6 




exp 


ik x /2 - \k y c y )W ■ W H r m r 




A 5 , 6 




exp 


ik x /2 - \k y c y )W ■ W H r m , 




A 6 ,6 




exp 


lk x /2 - ikyCy)W ■ W H , 




others 




o, 







(4.1) 



where 



= -1, 



= y/3/2, W = exp[-2(J x 



(3.11) 2J A2 + 2J B2 )/(k B T)}, W H = exp[-2H/(3k B T)], W a = 



exp[4JA2/(fc B T)], W b = exp[4J B2 /(fcBT 1 )], 
exp(i7r/6) and r m = exp(— iw/6). 

Then, the /^-function defined by (|3.3|) is 



D(k x , k y ) — M + c\ cosh(fc a; ) + 
Ci cosh(/ca;/2 — Cyky) + c\ cosh(fc :c /2 + c y k y ) + 
c 2 cosh(2/cj;) + c 2 cosh(fca; — 2c y k y ) + 
C2 cosh(/c a; + 2c y ky) + C3 cosh(2c y ky) + 
C3 cosh(3/ca;/2 — Cyk y ) + C3 cosh(3/ca;/2 + c y k y ) 

51 sinh(fc x ) + s 2 sinh(2fca;) — 

52 smh(k x — 2c y ky) + S4 sinh(fc K /2 — c y k y ) + 
S4 sinh(/c a; /2 + c y k y ) — S2 smh(k x + 2c y k y ), 



(4.2) 
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where 



M = 1 + 3W 2 + AW 6 - W 3 /W% - W 3 W\ 
12W e W a + 9W 6 W 2 + AW e W 3 + 
{W*W 3 )/W% + W 3 W^W 3 - 6W 6 W 4 
W 6 W® - l2W 6 W b + 30W 6 W a W b + 
{3W 3 W a W b )/W 3 H + 3W 3 WfjW a W b - 
l8W e W 2 W b - 6W 6 W 3 W b + 
6W 6 W 4 W b + 9W 6 W b - 



l8W 6 W a W b 2 

9w 6 w 2 wg -\ 

(W 3 W 3 )/W 3 
6W 6 W a W b 3 
6W 6 W b +6W 6 W a W l 



f 3w 4 wlwl - 

WV G W b 3 + 

W 3 W%W b 3 
2W & WlW b 3 - 



W e W, 6 



b ' 



Cl 



w 2 /w% - w^/wlj - 

W/W H - WW H + W 2 W 2 H - 
W 4 W^ + {W 4 W a )/W^ + W 4 W^W a 
(W 4 W 2 )/W 2 H 
{W 4 W 3 )/W 2 H 
{W 4 W b )/W 2 - 



- w 4 w 2 H wl 

- w 4 w 2 H wl 

W 4 W 2 H W b - 



(W 2 W a W b )/W 2 - (W 4 W a W b )/W 2 + 
(2W 3 W a W b )/W H - (2W 5 W a W b )/W H 
2W 3 W H W a W b - 2W 5 W H W a W b - 
W 2 W 2 H W a W b - 

W 4 W 2 H W a W b + (3W 5 W 2 W b )/W H + 
3W 5 W H W 2 W b - (W 5 W 4 W b )/W H - 

w 5 w H w 4 w b + (w 4 wl)jwl + 
w 4 w 2 H wl + 

{ZW^WIW^/Wh 

{w 4 w^)/wl - w 4 w 2 H wl - 

(W s W a Wf)/W H - W 5 W H W a W^ 



3W^W H W a Wl - 
- 3W 5 W H W% W b 2 



W H + WW% - W 3 Wji 



si = ~[W(-1 + W H )(l + W H )(W 

w 3 

W 3 W a + W 3 W 2 H W a 
W 3 W^W 2 - W 3 W 3 
W 3 W b + W 3 W^W b - 
ww a w b - W 3 W a W b - 
2W 2 W H W a W b + 2W 4 W H W a W b 



W 3 W 2 + 

w 3 w 2 H wl 



(4.3) 



(4.4) 



c 2 = [W 3 (l + W 2 H ){-1 + W a )(-1 + W b )]/W H , (4.5) 

c 3 = 2W 4 (-1 + W a ){-1 + W b )(-l + 

W a + W b + W a W b ), (4.6) 



ww%w a w b - W 3 W 2 H W a W b - 
3W 4 W H W 2 W b + W 4 W H W 4 W b 



w 3 wl + w 3 w 2 H wl 



W/ 4 W H WaWl 



?>W 4 WhWIwI 



S2 



Si 



w 3 w% - w 3 w 2 H wl 
w 4 w H w a wt)\lw 2 H , 



-[W 3 {-l + W H ){l + W H ){-l 
+W a ){-1 + W b )]/W H , 



(4.7) 



(4.8) 



[W{- 

w 3 - 

W 3 W a - 

W 3 W 2 

W 3 W 3 



■1 + W H )(1+W H )(W 
W H + WW% - W 3 W 



■ H 



W 3 WfjW a - 

- w 3 w 2 H wl 

■ W 3 WhW 3 



W 3 W b + W 3 W 2 H W b - 

WWaWh - W 3 W a W b - 

2W 2 W H W a W b + 2W 4 W H W a W b 
WWjjW a W b - W 3 WjjW a W b - 
3W 4 W H W 2 W b + W 4 W H W 4 W b - 



W 3 Wi 



w 3 w 2 H wl 



2>W 4 W H W a Wl 

w 3 wf 



m 4 w H w 2 a wl 



w 4 w H w a wt)\lw 2 H . 



(4.9) 



We substitute the D-function into (3.4) and solve them 
with respect to {uj x ,uj y ) as a function of 9. Substituting 
the solution (ui x (8), oj y (9)) into (3.5)— (3.8), we obtain the 
interface tension, 2D island shape and the interface stiff- 
ness. 

Note that the Z?-function has the mirror symmetry 
with respect to 

Therefore, the island shape has the mirror symmetry 
with respect to fc^-a xis. That is, uj y {Q) — and u y (-K) = 
are the solutions of (3.4). At the orientation correspond- 
ing to 6 = or 9 — ir, the form (4.2) reduces to 

D(k x , 0) = M + c 3 + ci cosh(fc ;E ) + 
2ci cosh(/c 2; /2) + +c 2 cosh(2fca;) + 
2c2 cosh(/c :E ) + 2c3 cosh(3fc x /2) + 
si sinh(fc K ) + S2 sinh(2fc x ) — 
2s 2 sinh(fe x ) + 2s 4 sinh(fc 2 ,/2). (4.10) 

From the solution of D(k x ,0) — (u>y(0) — u y (ir) = 0), 
we o btain cosh(wa;(0)/2) and cosh.(aj x (w)/2). Then, from 
,5[), step tension becomes 



7 (0) = 2fc B Tcosh- 1 (^(0)/2), 
7(tt) = 2/c B 7 1 cosh~ 1 (u; K (7r)/2). 



(4.11) 



In the T — > limit, step tensions (step free energy per 
lattice constant) becomes, 
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7 (0) = Min[2 Jx + 4J 2 + -H, 



B. The case of H = 



2(2Ji+ 47 a --#)], 
7 (tt) = Min[2 Ji + 4J 2 - | fl, 



2(2 Ji + 4J 2 



(4.12) 



where Min[a, b] denotes the smaller one in {a, b}. 
From (3.7), the step stiffness becomes, 



7(0) = AfeT[|ci sinhK.(0)/2) + 
(ci + 2c 2 )smh(w x (0)) + 
3c 3 sinh(3w x (0)/2) + 2c 2 sinh(2a; x (0)) + 
s 4 cosh(a; x (0)/2) + (s x + 2s 2 ) cosh(w 2 .(0)) + 
2s 2 cosh(2 W:E (0))|]/[2c 2 (2c 3 + 
ci cosh(a; x (0)/2) + 4c 2 cosh(u> x (0)) + 
c 3 cosh(3u; x (0)/2) + s 4 sinh(a; x (0)/2) - 
4s 2 smhK(0)))], (4.13) 

7(tt) = fc B T[| Cl sinh(w x (7r)/2) + 
(ci + 2c 2 ) sinh(u; x (7r)) + 
3c3 sinh(3<J x (7r)/2) + 2c 2 sinh(2u> x (7r)) + 
s 4 cosh(w x (7r)/2) + (si + 2s 2 ) cosh(w x (7r)) + 
2s 2 cosh(2a; x (7r))|]/[2c^(2c3 + 
ci cosh(u; x (7r)/2) + 4c 2 cosh(a; x (7r)) + 
c 3 cosh(3wic (7r) /2) + s 4 sinh(a) x (n) /2) — 
4s 2 sinhK(7r)))]. (4.14) 

In Fig. |^, we show an example of equilibrium island 
shape and a polar graph of step stiffness at J\ = 165meV, 
J 2 /J x = -0.1, lattice constant= 3.84A and H/J x = 0.31. 
We also show the temperature dependence of step ten- 
sion, step stiffness, the coefficient of step interaction 
( |3.1l[) where g = B/a\, a h = 3.14A and g = 0. 

In the absence of nnn interactions, W a and Wb reduce 
to unity. The D-function (4.2), then, becomes 

D(k x , k y ) = M + ci cosh(fe x ) + 
2ci cosh(fc x /2 — Cyky) + ci cosh(fe x /2 + c y k y ) + 
S\ sinh(k x ) + s 4 sinh(fc x /2 — c y k y ) + 
s 4 sinh(fc 2; /2 + c y k y ), (4-15) 

M = 1 + 3W 2 + 3W 4 + W 6 + AW Z {W% + 

ci = -(1 - WfW{\ + Wf{l + W%)/W B , 
si = -s 4 

= (1 - W) 2 W(1 + VT) 2 (1 + Wh)(1 - 
W H )/W H , 

c 2 = c 3 = s 2 = 0, (4.16) 



At 77 = 0, si, s 2 and s 4 in (4.2) become zero, since 
Wh reduces to unity. The £>-function, then, has the 
symmetry D(k x ,k y ) — D(±k x ,±k y ). The island shape 
has the mirror symmetry with respect to the fc y -axis too; 
ui x (ir/2) = u> x (3ir/2) = becomes the solution of (|3.4|). 
The equation 

£>(0, k y ) = M + ci + c 2 + 2(ci + c 3 ) cosh(c y fc y ) + 
(2c 2 + c 3 ) cosh(2c y /c y ) 
- 0, (4.17) 

is solved, in terms of cosh(c y ujy) {k y — iui y ), as 

ci + c 3 



cosh(c y Wj / ) 



2(2c 2 + eg) 



(ci + c 3 ) 2 — ci + c 2 + c 3 - M 



V 4(2c 2 + c 3 ) 2 2(2c 2 + ca) 

= z. (4.18) 

By use of this solution, we obtain step tension as 

7(tt/2) = 7 (37t/2) 

■ cosh _1 (z), 



(4.19) 



and the step stiffness as 
7 (tt/2) = 7 (3tt/2) 



-4k B Tc y yjz 2 - l\ Cl + c 3 + 2(2c 2 + c 3 )z| 
2ci + 4c 2 + (ci + 9c 3 )z + 8c 2 z 2 



(4.20) 



In the case of Ja2 = Jb2 = 0, the Z?-function reduces 
to that of the exact solution for the nn honeycomb lattice 
system [f22f . 

That is, 

D{k x , k y ) = M + Cj cosh(fc x ) + c 4 cosh(fc x /2 — c^fe^) + 
Ci cosh(fc,./2 + Cyky), 

M = (1 + VK) 2 (1 - 2W + 6W 2 - 

2W 3 + W 4 ), 
ci = -2(1 - W0 2 W(1 + W) 2 , 
c 2 = c 3 = si = s 2 = s 4 = 0. (4.21) 

At 8 = 0, we obtain an explicit form of the solutions as 

LOy(0) = 0, 

cosh(cj :c (0)/2) = ~ ^3 - 2M/ci - i. (4.22) 
Also, at 6 — 7r/2, we have 



which agrees with the £>-function given in Rcf. 1 1 3 1 



cj x (tt/2) = 0, 
cosh(c^(V2)) = \(-M/ Cl - 1). 



(4.23) 
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Hence, the interface tensions become 

7 (0) = 2fc B Tln(m. y ), 
1 1 



= -77 + ^2 + - A /( Z2 -3)(z 2 + l), 



Z 2 



(i + w)Vi-w + w 2 



W\l - w\ 



(4.24) 



7 (tt/2) = 1/cy ■ fc B Tln[(-2 + 1/W + W)/2]. (4.25) 

Due to the Wulff's theorem, 7(0) and 7(7r/2) give the 
linear size of the island shape along x- and y-direction. 
The explicit form of interface stiffness becomes 



7(0) 



fc B T - 2z 2 + z 2 



c 2 y 2(z 2 - 1) 
7(tt/2) = Ak B Tc y 

{1 + W 2 )\l - AW + W 2 \ 



1 + AW ~ 6W 2 + AW 3 + W 4 ' 



(4.26) 



(4.27) 



2 l~W 2 
73(0) = ^eTlnh^H 

73 (7r/2) = 2fc B T x 



cosh 1 ( — ~ + r 4 / 3 



1 + 3IU 4 



2 2y VK 2 (1 - W)(l + W) 
The step stiffness is written as 

2v / 3fc B T(l + V^ 2 )|l-3W^ 2 | 



73(0) = 



1 + 6W 2 - 3W 4 



(4.30) 



(4.31) 



(4.32) 



73(tt/2) = — — x 



V(l + 3W 2 )(1 + 3VF 4 - 2Wz x + 2W 3 Zl ) 



W(zi - W)\l - W 2 \ 



z\ = 



1 + 3W 2 
1 - W 2 ' 



(4.33) 



V. APPLICATION TO SI(lll) SURFACE 



C. Triangular lattice 

In the limit of J\ — > 0, the system becomes two inde- 
pendent triangular lattice gases. In this case, the random 
walk treatment on the honeycomb lattice may not be 
good enough. In fact, (4.2) in this limit gives D-function 
slightly different from the known exact one on the trian- 
gular lattice. Hence, we need a separate study to treat 
this case. 

The D-function of the triangular lattice becomes 

D 3 (k x ,k y ) — M — c\ cosh(fc a ) — c 2 cosh(v / 3fc a; /2 — k y /2) — 
c 3 cosh(V3k x /2 + k y /2), (4.28) 



where 



M = 1 + W\WI + WlWl -I- W%W? 



3 > 



ci = 2(1 - + W 1 )W 2 W 3 , 

c 2 = 2W 1 W 2 (1 - W 3 )(l + W 3 ), 
c 3 = 2W 1 (l-W 2 ){l + W 2 )W 3 , 



(4.29) 



Wi = exp[-2Ji/(fc B T)], W 2 = exp[-2J 2 /(fc B T)], and 
W 3 — exp[— 2 J 3 /(fc B T)]. The equations agree with the 
known exact one p2fl . Note that the form of D 3 (k x , k y ) 
(4.27) is essentially the same as D(k x ,k y ) of nn honey- 
comb lattice (4.20). Therefo re, t he island shape of trian- 
gular lattice obtained from (3.4) is the same as the one 
of honeycomb lattice. The difference is the temperature 
dependence of coefficients. When J\ = J 2 = J 3 = J, 
i.e., Wi = W 2 — W 3 — W, ZJ-function has a symmetry 
of D 3 (k x ,k y ) = D 3 (±k x ,±k y ). Hence, we have explicit 
forms of 7 and 7 for special orientations. Therefore, we 
have 



A. The 7x7 reconstructed surface 

At temperatures lower than the 7x7^ lxl transition 
temperature (~1130 K), Si(lll) surface forms 7x 7 recon- 
structed structure called DAS (dimer adatom stacking- 
fault) structure |23j . The unit cell of the DAS structure 
is divided into the faulted half (FH) and the unfaulted 
half (UH) pi} |, each of which forms a triangular lattice. 
It has been observed that the step structure is well de- 
scribed by the combination of the FH unit and UH unit 
p5| . We consider, therefore, a pair of triangular sub- 
lattices, where the one represents the FH lattice system, 
and the other represents the UH one. Consequently, the 
system becomes equivalent to a stoichiometrically binary 
lattice gas on a honeycomb lattice with nn interactions, 
where the inequivalent sites of the lattice gas model are 
coarse-grained representations of these two halves of the 
7x7 unit cell. Therefore, we can use (4.15) and (4.16) 
to calculated step quantities. In Fig. ||, we regard closed 
circles as FH units, and open circles as UH units. We set 
lattice constant = 3.84 x 7 A, and step height = 3.14A. 
We introduce the step running direction angle 9 so that 
a straight step with 9 = corresponds to (112) step (Fig. 

In spite of the extensive experimental studies on 
Si(lll)7x7 structure, the values of kink energy and step 
tension have not been settled yet. As one trial, we adopt 
the result of Eaglesham et al. |2(| where step tension 7 at 
700 °C was obtained from the equilibrium crystal shape 
(ECS) of Si: 7111 = 5.7 x 10~ n J/m for (111) surface 
and 7100 = 1.0 x 10~ n J/m for (100) surface. Due to 
the 2x1 reconstruction, this value of 7100 corresponds 
to the mean value of the S^-step tension and the S's-step 
tension, and is consistent with the one calculated in our 
previous papers HJlT 
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We choose kink energy so that the calculated mean 
value of step tension for (211) and (211) at 700 °C re- 
produces the above mentioned value 5.7 x 10 -11 J/m 
(= 36 meV/A). In addition, the experimental observa- 
tion of the island shape (and also the shape of spiral 
step) gives further information on the kink energy, due 
to the Wulff's theorem for 2D ECS. Let h n be the dis- 
tance between the center (= Wulff point) of the ECS 
(island shape, in our case) and the tangential line of the 
ECS at a position on the ECS curve, where n is the 
interface normal vector at the position. The Wulff's the- 
orem states that the ratio hn/~fn {in- interface ten- 
sion, or step tension in our case) is n-independent, lead- 
ing to a relation jn/jrv — hn/hn' for arbitrary n 
and n' . From the photographs of the experimental ob- 
servation p7j , p8|] , we have fcm/^il = 1-2. This ra- 
tio gives the ratio between the step free energies cor- 
responding to these directions. At such low temperature 
where the observation was made, these step free ener- 
gies are well approximated by the step formation ener- 
gies (2J + 2H/3) and (2J - 2H/3) (see (4.12)), giving 
us /1211//1211 = I- 2 = (2J + 2H/3)/(2J - 2H/3) which 
amounts to H/J = 0.31. 

We then set J = 0.475eV and H = 0.15eV. The kink 
energy becomes 1.05eV= 2J + 2H/3 for (211) step and 
0.85eV= 2J - 2H/3 for (211) step, which are smaller 
than but is in the same order of magnitude of the ones 
reported in Refs. |^@. 

The difference in the on-site energy between the UH 
and FH, is then E FH - E UH = 4H = 0.59eV. In the first 
principles study of Meade and Vanderbilt |3C[ ] , surface 
energies of the Si(lll) surface for various structures are 
calculated: For example, 1.24 eV/1 x 1 for the 2x2- 
adatom structure, and 1.27 eV/1 x 1 for the 2 x 2-adatom 
(faulted) structure. From these values, we can estimate 
E FH -E UH to be (1.27- 1.24) x 24 - 0.7eV which is in 
reasonable agreement with our value 0.59eV. 

In Fig. we show equilibrium island shape at 400 
°C and 850 °C, and the temperature dependence of 
step tension, step st iffnes s and step interaction coeffi- 
cient g = B/a^ (see ( 3.11 )). The step tension is almost 
constant below 1130K, because the temperature is very 
low as compared with the lattice-gas melting tempera- 
ture of the model (~ 8300K). On the other hand, the step 
stiffness strongly depends on temperature in the same re- 
gion. The step stiffness for (211) step increases rapidly as 
temperature decreases. While, the step stiffness of (011) 
step becomes smaller and smaller and converges to zero 
at zero temperature. 

Note that a similar analysis can be made for n x n DAS 
structure. From the photographs of small island of 5 x 5 
structure (28), we find that the island shape has a six- 
fold rotational symmetry in contrast to the case of the 
7x7 structure which has a three-fold rotational symme- 
try. Recall that the six-fold rotational symmetry appear 
only when H = (see §4.2). Therefore, the energy differ- 
ence between FH and UH units in the case of 5 x 5 DAS 
structure is very small if exists. We stress here that, also 



for other structures, observation of the anisotropy of the 
equilibrium island shape will be useful in determining the 
energy difference between FH and UH units. 



B. The lxl high temperature surface 

The high-temperature Si(lll) surface at about 900 °C 
has the structure of the lxl surface together with disor- 
dered adatoms with concentration of 0.25 ]3l|-|33]|. Fur- 
ther, Khomoto and Ichimiya |3lJ reported that the num- 
ber ratio of adatoms sitting on T4 site and H3 site is 
4:1. Although the adatoms are considered to be in a 
disordered phase, the broad \/3 x \/3 peaks appear in 
diffraction observations (3^j3l[] which suggests existence 
of the short-range order corresponding to formation of 
the hard-hexagon units Q . 

For the high temperature surface of Si(lll), experi- 
mental measurement of the step tension and the step 
stiffness has been a subject of active study |^l|,^-^9). 
The experimental values are, however, not settled yet. 
We make, therefore, several trial calculations for possible 
cases. In all the cases, we choose the microscopic cou- 
pling constants so that the calculated step stiffness at 
900°C reproduces the value presented by Bartelt et al. 



1. Case 1: Adatom in disordered phase without short-range 

order 

We use the honeycomb lattice gas system of (4.2), with 
J\ and J2 being regarded as effective coupling constants. 

In Fig. ^| we regard the filled circles as atoms of the top 
layer and the open circles as those of the second layer (the 
lattice constant = 3.84A, the step height — 3.14A). 

We introduce the step running direction angle 9 so that 
a straight step with 9 — corresponds to (112) step. The 
effect of dangling bonds normal to (111) plane are taken 
into account by setting H/Ji = 1. 

We calculate equilibrium island shape, step stiffness 
and step interaction coefficient, which we show in Fig. |(]. 
Here, assuming that J 2 is small, we have set J2/J1 = 0.2 
and J\ — 60meV. The kink energy becomes 176meV for 
(211) step and 256meV for (211) step. As is seen in 
the island shape and the polar graph of step stiffness at 
900°C, there remains anisotropy in the step stiffness in 
spite of the circular island shape. The most significant 
characteristic of this results is the asymmetry between 
the orientations (211) and (211). In contrast to the 7x7 
structure, the stiffness as a function of the step orienta- 
tion takes its maximum at (211). 
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2. Case 2: The x ^3 short -range ordered phase 

In the case of \/3 x \/3-ordered phase of adatoms, we 
calculate step quantities by using the triangular lattice 
gas model (Fig. |7|a). As has been pointed out in Sec. 4, 
the system has six-fold rotational symmetry. We set the 
lattice constant to be 3.84 x v^A, and the step height 
to be 3.14A. We introduce the step running direction 
angle 9 so that a straight step with 6 — corresponds 
to (112) step. The step tension and the step stiffness 
are calculated exactly from (4.30)-(4.33). The effective 
coupling constant is obtained to be J = 62meV (kink 
energy = 248meV). We show the calculated results in 
Fig. ||. The step stiffness takes its maximum at the 
orientation {10l}. 



3. Case 3: The 2x2 short-range ordered phase 



with the numerical renormalization-group method; 
details will be discussed elsewhere. 
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Calculation of step quantities can be done in the same 
fashion as in the case 2. We set the lattice constant as 
3.84 x 2A, and the step height as 3.14A (Fig. 0b). We 
introduce the step running direction angle 9 so that a 
straight step with 9 = corresponds to (Oil) step. The 
step stiffness takes its maximum at {211}. The effective 
coupling constant is obtained to be J = 67meV (kink 
energy = 268meV). We show the calculated results in 
Fig- H 

VI. SUMMARY 

We have considered the honeycomb lattice Ising sys- 
tem in a staggered field with both nearest-neighbor (nn) 
and next-nearest-neighbor (nnn) interactions, to calcu- 
late interface tension, interface stiffness, island shape by 
the imaginary path-weight (IPW) method. 

We have applied the calculated results to Si(lll) 7x7- 
reconstructed surfaces and the high-temperature Si(lll) 
lxl surface. We have made estimation on the mi- 
croscopic coupling constants from existing experimental 
data, and have drawn equilibrium island shape, step ten- 
sion, step stiffness and the coefficient of step interaction, 
with their temperature dependence. Our analysis made 
in the present paper will be helpful in determining precise 
value of the kink energy from experimental observation. 

Our lattice-gas treatment made in the present paper 
corresponds to the two-level approximation for the sur- 
face fluctuation. Fortunately, the temperature-range of 
our concern in the present study is very low, the two-level 
approximation is expected to be fairly reliable. On the 
other hand, at higher temperatures, near the roughen- 
ing transition temperature, we should consider multilevel 
fluctuation of the surface. Even in such cases, we have an 
efficient method, namely, the temperature-rescaled Ising- 
model approach, Jl5[ where the IPW method is combined 
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FIG. 1. Examples of an interface of square lattice Ising 
model made by fixing the boundary spins. 



FIG. 2. Examples of an interface configuration. A-atom 
and B-atom are indicated by filled circle and open circle, re- 
spectively. Thick line represents an interface. 



FIG. 3. An example of calculation by the use of the 
D- function of (4.2). (a) The island shape at 900° C, (b) a 
polar graph of step stiffness at 900° C, (c) temperature de- 
pendence of step tension, (d) temperature dependence of step 
stiffness, and (e) temperature dependence of g — Bja\. We 
have set Ji = 165meV, J2 = — 16.5mev, and H — 165meV: 
Kink energies are 88meV for (211) step and 176meV for (211) 
step. In (c)-(e), thick lines correspond to(211) step, thin lines 
to (211) step and broken lines to {101} step. 



FIG. 4. Examples of a step edge on 7 x 7 reconstructed 
Si(lll) surface. "U" denotes unfaulted half unit, and "F" 
denotes faulted half unit. 



FIG. 5. Calculation for 7 x 7 reconstructed surface by 
the use of the D- function of (4.15). (a) The island shape at 
850°C, (b) the island shape at 400°C, (c) temperature de- 
pendence of step tension, (d) temperature dependence of step 
stiffness and (e) temperature dependence of g = -B/af . We 
have set J = 0.475eV, 4H = 0.59eV: Kink energy = 1.05eV 
for (211) step and 0.85eV for (211) step. In (c)-(c), thick lines 
corresponds to (211) step, thin lines to (211)step and broken 
lines to {101} step. 



FIG. 6. Calculation for 1 x 1 surface (Case 1) by the 
use of D-function of (4.2). (a) The island shape at 900°C, 
(b) a polar graph of step stiffness at 900° C, (c) temperature 
dependence of step tension, (d) temperature dependence of 
step stiffness and (e) temperature dependence of g — B/a\ . 
We have set Ji = 60meV, J 2 = 12mev, and H = 60meV: Kink 
energy = 176meV for (211) step and 256meV for (211) step. 
In (c)-(e), thick lines corresponds to (211) step, thin lines to 
(211)step and broken lines to {101} step. Open squares: Ref. 
[37]. Open circle: Ref. [35]. 
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FIG. 7. Examples of an interface configuration for two 
cases of adatom orderings: (a) \/3 x \/3 and (b) 2x2. Thick 
line represents the interface. A-atom and B-atom are indi- 
cated by filled circle and open circle, respectively. Adatoms 
are indicated by shaded large circles. Broken lines denote 
boundaries between the hexagons. 



FIG. 8. Case 2. Calculation for ^3 x ^3 adatom order- 
ing (Case 2) by the use of D- function of (4.30)-(4.33). (a) 
The island shape at 900° C, (b) a polar graph of step stiffness 
at 900°C, (c) temperature dependence of step tension, (d) 
temperature dependence of step stiffness, (e) temperature de- 
pendence of g = B/a\. We have set J = 62meV: Kink energy 
= 248meV for {101} step. In (c)-(e), thick lines corresponds 
to {211} step and thin lines to {101} step. Open squares: 
Ref. [37]. The open circle: Ref. [35]. 



FIG. 9. Case 3. Calculation for 2x2 adatom ordering 
(Case 3) by the use of D-function of (4.30)-(4.33). (a) The 
island shape at 900° C, (b) a polar graph of step stiffness at 
900° C, (c) temperature dependence of step tension, (d) tem- 
perature dependence of step stiffness and (e) temperature de- 
pendence of g = B /a\. We have set J = 67meV: Kink energy 
= 268meV for 211 step. In (c)-(e), thick lines correspond to 
{211} step and thin lines to {101} step. Open squares: Ref. 
[37]. The open circle: Ref. [35]. 
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Figure 9 N. Akutsu and Y. Akutsu 



